{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The proof of the [Perron-Frobenius](https://en.wikipedia.org/wiki/Perron%E2%80%93Frobenius_theorem) theorem can seem very abstract, but if you play with some examples it is easier to understand.\n",
    "This notebook presents the proof with computational examples.<br>\n",
    "\n",
    "Step 4 below uses JuMP to turn Perron-Frobenius into a computational algorithm.\n",
    "<br>\n",
    "\n",
    "There are a few variations on the theorem some with more and some with less information\n",
    "but the basic version says that if A has all positive entries, then the maximum\n",
    "absolute eigenvalue is real and positive and there is a corresponding real positive eigenvector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step #1. Assume all(x.>0) and all(y.>0) and define τ as the minimum of y./x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "τ (generic function with 1 method)"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Define τ(y,x) on vectors\n",
    "\n",
    "τ(y::Vector, x::Vector) = minimum(y./x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that for 0 ≤ t ≤ τ(y,x) we have all(y .≥ t*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Example\n",
    "y = [10,5,6,9]\n",
    "x = [1,2,3,4]\n",
    "τ(y,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(true, true, false)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "all(y.≥2x), all(y.≥1.99x),all(y.≥2.01x)  # check these by hand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step #2. If all(A.>0) and  all(y.≥0) and y is not the zero vector then all(A*y.>0)  (strictly greater)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 0.2\n",
       " 0.5\n",
       " 0.8"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Example\n",
    "A= [ 1 2 3;4 5 6; 7 8 9]\n",
    "y = [0, .1, .0]\n",
    "A * y # any one positive entry multiplies an entire positive column of A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step #3: <br> τ(Ax,x)=τ(A²x,Ax) if x in an eigenvector with all(x.≥0). <br>\n",
    "τ(Ax,x) < τ(A²x,Ax) if x is not an eigenvector.\n",
    "\n",
    "<br>\n",
    "Proof:  If x is an eigenvector, then τ(Ax,x)=τ(A²x,Ax)= the corresponding eigenvalue. <br>\n",
    "If x is not an eigenvector, then letting y\n",
    "= Ax - τ(Ax,x) *x, then  all(y.≥0)  and y is not the 0 vector. <br>\n",
    "From Step 2, all(A*y.>0) or equivalently all(A²x .>  τ(Ax,x) *Ax) from which we see\n",
    "τ(A²x,Ax) > τ(Ax,x)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7-element Array{Float64,1}:\n",
       " 1.34884\n",
       " 2.40402\n",
       " 2.68214\n",
       " 2.75781\n",
       " 2.77293\n",
       " 2.77552\n",
       " 2.77666"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# An example\n",
    "n = 6\n",
    "A = rand(n,n)\n",
    "x = rand(n)\n",
    "[τ(A^k*x, A^(k-1)*x) for k=1:7] # This sequence will be increasing, but to an eig limit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step #4.  Let tmax be the maximum of τ(Ax,x) for all non-zero x. We will prove mathematically that x is a positive eigenvector and τ(Ax,x) is the eigenvalue. Before we do it mathematically, let's see it computationally:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way to form this maximum problem is write this as a constrained optimization:\n",
    "\n",
    "$\\max t$ subject to  <br> <br>\n",
    "$x_i \\ge 0 $ <br>\n",
    "$y=Ax$ <br>\n",
    "$y[i]/x[i] \\ge t$ <br>\n",
    "$sum(x)=1$\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use the highly popular Julia [Jump Package](https://github.com/JuliaOpt/JuMP.jl) created at MIT (though not in math!), and used widely for operations research and in business schools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pkg.add(\"JuMP\")\n",
    "using JuMP\n",
    "# Pkg.add(\"Ipopt\")  (On my mac, this worked with 0.6.2 but not 0.6.0)\n",
    "using Ipopt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7×7 Array{Float64,2}:\n",
       " 0.971603   0.325743  0.863038  0.0234046  0.962918   0.496618  0.799348\n",
       " 0.62474    0.140906  0.448296  0.505187   0.0646877  0.149136  0.205624\n",
       " 0.448767   0.58146   0.47302   0.443701   0.303789   0.114217  0.892493\n",
       " 0.808785   0.588347  0.839119  0.883789   0.920193   0.515088  0.22442 \n",
       " 0.0089511  0.242133  0.783681  0.420531   0.965035   0.544011  0.334241\n",
       " 0.300799   0.990369  0.401669  0.427284   0.207415   0.309122  0.329326\n",
       " 0.0314547  0.723179  0.476076  0.445037   0.249261   0.404243  0.502455"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(7,7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.3686909584508244"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n=size(A,1)\n",
    "\n",
    "m = Model(solver=IpoptSolver(print_level=2))\n",
    "@variable(m, t);         @objective(m, Max, t)\n",
    "\n",
    "@variable(m, x[1:n]>=0); @constraint(m, sum(x)==1)\n",
    "@variable(m, y[1:n]);    @constraint(m, y .== A*x)\n",
    "\n",
    "@NLconstraint(m, [i=1:n], t <= y[i]/x[i])  # nonlinear constraint\n",
    "\n",
    "status = solve(m)\n",
    "x = getvalue.(x)\n",
    "t = getobjectivevalue(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.4817758919083086e-6"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(A*x-t*x) #  demonstrate we have found an eigenpair through optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 5: Demonstrate that if x above were not an eigenvector, then the t could not have been the solution to the optimum problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we saw in step 3, if x had not been an eigenvector, then τ(Ax,x) < τ(A²x,Ax), so τ(Ax,x) was not the maximum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Step 6:  Any complex eigenvector, eigenvalue pair has absolute eigenvalue <= tmax:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If Ax = λx then all( A*abs.(x) .≥ abs(λ)*abs.(x)) by the triangle inequality.  Thus abs(λ) <= tmax."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×5 Array{Float64,2}:\n",
       " 0.996711  0.656579  0.453247   0.61344   0.50697  \n",
       " 0.591166  0.616613  0.987583   0.246784  0.442663 \n",
       " 0.949881  0.454748  0.831274   0.708647  0.458239 \n",
       " 0.069995  0.108182  0.0296905  0.434673  0.0322304\n",
       " 0.105186  0.918176  0.831151   0.126704  0.0709903"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(5,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Complex{Float64},1}:\n",
       "   2.58617+0.0im    \n",
       "  0.125586+0.34277im\n",
       "  0.125586-0.34277im\n",
       "  0.351225+0.0im    \n",
       " -0.238306+0.0im    "
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigvals(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1255862966495158 + 0.3427698069874712im"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Λ,X=eig(A);x=X[:,2];λ=Λ[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.355107029738416e-16"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(A*x-λ*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6784932085048402"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "τ(A*abs.(x),abs.(x)) - abs(λ) # This is non-negative"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.2",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
